#include <cmath>
#include <gsl/gsl_randist.h>
#include "networkplanning.h"
#include "channel.h"
#include "simulator.h"

extern const gsl_rng_type * T;
extern gsl_rng *r;

double
CalculatePathLoss(int x0, int y0, int x, int y, bool LOS)
{
	// Calculate path loss between different block
	// Reference: 802.16m p.40 3.2.3.3
	double remain;
	bool flag;
	double p1 , p2 , R , r0 = -1.0;
	double p;

	double distance = sqrt( pow( double( PRECISION * ( x0 - x ) ) , 2 )
			              + pow( double( PRECISION * ( y0 - y ) ) , 2 ) );
	if( distance == 0 )
	{
		return 0;
	}
	if( distance <= 15 || LOS )
	{
		flag = true;
	}
	else
	{
		remain = 1 - pow( 1 - pow( 1.56 - 0.48 * log10( distance ), 3 ) ,
				          1 / 3.0 );
		if( gsl_rng_uniform (r) <= remain )
		{
			flag = true;
		}
		else
		{
			flag = false;
		}
	}
	if( flag )  //line of sight(LOS)
	{
		p = 32.4418 + 20 * log10 ( FREQUENCY )
		    + 20 * log10( distance ) + .0174 * distance
		    + 20 * log10( ( 0.013 * distance ) > 1 ? ( 0.013 * distance / FREQUENCY ) : 1 );
		return p;
	}
	else        //non-line of sight(NLOS)
	{
		R = PRECISION * ( abs( x - x0 ) + abs( y - y0 ) );
		p1 = 24 + 45 * log10( distance + 20 );

	    //int segments = round( distance / (sqrt(2) * 50) );
		//Defined in the specification, but without practical use in our calculation,

		// a practical solution, not necessarily right
		double k = 1, d = 0,
			   dx = PRECISION * abs( x - x0 ), dy = PRECISION * abs( y - y0 );
		dx -= 50; dy -=50;    //50m is the segment length
		while (dx > 0 && dy > 0)
		{
			r0 = 50.0;
			k += d;
			d += k * 50;
			k += d;
			d += k * 50;
			dx -= 50; dy -=50;
		}
		k += d;
		if (dx < 0 && dy < 0)
		{
			d += k*(50+dx);
			k += d;
			d += k*(50+dy);
			if (r0 < 0)   //Caution: r0 is double, cannot use "r0==50"
			{
				r0 = (dx > dy) ? (50+dx) : (50 + dy);
			}
		}
		else if (dy >= 0)
		{
			d += k*50;
			k += d;	d += k*(50+dx);
			k += d;	d += k*dy;
			if (r0 < 0)
			{
				r0 = 50;
			}
		}
		else
		{
			d += k*50;
			k += d;	d += k*(50+dy);
			k += d;	d += k*dx;
			if (r0 < 0)
			{
				r0 = 50;
			}
		}
		double rbp = (76.67*FREQUENCY  < r0 ? 76.67*FREQUENCY : r0 );

		p2 = 32.4418 + 20 * log10 ( FREQUENCY ) + 20 * log10( d )
			 + 20 * log10( ( R / rbp ) > 1 ?  (R / rbp) : 1 )
			 + 20 * log10( R ) + .0174 * R;

		p = p1 < p2 ? p1 : p2;
		return p;
	}
}
